# -*- coding: utf-8 -*-

"""
@File    : 03_读取栅格像元坐标.py
@Author  : fungis@163.com
@notice  :
"""

from osgeo import gdal


def write_txt_by_list(path, list=None):
    '''

    :param path: 储存list的位置
    :param list: list数据
    :return: None/relist 当仅有path参数输入时为读取模式将txt读取为list
             当path参数和list都有输入时为保存模式将list保存为txt
    '''
    if list != None:
        file = open(path, 'w')
        file.write(str(list))
        file.close()
        return None
    else:
        file = open(path, 'r')
        rdlist = eval(file.read())
        file.close()
        return rdlist


if __name__ == '__main__':
    filePath = r'./data-use/tif/NDVI_201405_bm.tif'  # 输入你的栅格数据
    dataset = gdal.Open(filePath)
    adfGeoTransform = dataset.GetGeoTransform()
    # 打开波段1（注意:用索引1，而不是0，来获取第一个波段）
    band = dataset.GetRasterBand(1)
    band_array = band.ReadAsArray()  # 一维像元矩阵
    # 左上角地理坐标
    print(adfGeoTransform[0])
    print(adfGeoTransform[3])

    nXSize = dataset.RasterXSize  # 列数
    nYSize = dataset.RasterYSize  # 行数

    location_list = []  # 用于存储每个像素的（X,Y）坐标
    cell_list = []  # 用于存储每个像素的（X,Y,Z）坐标
    for i in range(nYSize):
        location_row = []
        cell_row = []
        for j in range(nXSize):
            px = adfGeoTransform[0] + i * adfGeoTransform[1] + j * adfGeoTransform[2]
            py = adfGeoTransform[3] + i * adfGeoTransform[4] + j * adfGeoTransform[5]
            location_row.append([px, py])
            cell_col = [px, py, band_array[i][j]]
            cell_row.append(cell_col)
        location_list.append(location_row)
        cell_list.append(cell_row)
    print(len(location_list))
    write_txt_by_list(r'./results/test221012.03.txt', location_list)
    write_txt_by_list(r'./results/test221012.04.txt', cell_list)
